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ABSTRACT 

Observations suggest that star clusters often form in binaries or larger bound groups. Therefore, mergers between two clusters 
are likely to occur. If these clusters both harbor an intermediate-mass black hole (IMBH; ~ IO^Mq) in their center, they can 
become a strong source of gravitational waves when the black holes merge with each other. In order to understand the dynamical 
processes that operate in such a scenario, one has to study the evolution of the merger of two such young massive star clusters, 
and more specifically, their respective IMBHs. We employ the direct-summation NBODY4 numerical tool on special-purpose 
GRAPE6 hardware to simulate a merger of two stellar clusters each containing 63,000 particles and a central IMBH. This allows 
us to study accurately the orbital evolution of the colliding clusters and the embedded massive black holes. Within ~ 7 Myr the 
clusters have merged and the IMBHs constitute a hard binary. The final coalescence happens in ~ 10 8 yrs. The implication of our 
analysis is that intermediate-mass black holes merging as the result of coalescence of young dense clusters could provide a source 
for the Laser Interferometer Space Antenna (LISA) space-based gravitational wave detector mission. We find that interactions 
with stars increase the eccentricity of the IMBH binary to about 0.8. Although the binary later circularizes by emission of 
gravitational waves, the residual eccentricity can be detectable through its influence on the phase of the waves if the last few 
years of inspiral are observed. For proposed higher-frequency space-based missions such as the Big Bang Observer (BBO), 
whose first purpose is to search for an inflation-generated gravitational waves background in the 10 _1 — 1 Hz range, binary IMBH 
inspirals would be a foreground noise source. However, we find that the inspiral signals could be characterized accurately enough 
that they could be removed from the data stream and in the process provide us with detailed information about these astrophysical 
events. 

Subject headings: Black hole physics, gravitational waves, stellar dynamics, methods: N-body simulations 



1. INTRODUCTION 

There are at least two lines of evidence indicating that stel- 
lar clusters, containing 10 4 — 10 6 stars, can merge with each 
other relatively early in there evolution. First, high-resolution 
Hubble Space Telescope observations of the Antennae (Whit- 
more et al. l!999HZhang & Falllll999l) or Stephan's Quintet 
JGallagheretal.120011) reveal hundreds of young massive star 
clusters in the star forming regions. These clusters are clus- 
tered into larger complexes of a few 100 pc. Since they har- 
bor ~ 10 5 stars within a few parsecs and are older than 5 Myr, 
they are most likely bound clusters. These cluster complexes 
have been suggested as the progenitors of ultra-compact dwarf 
galaxies (UCDGs), as a result of t he amalgamation of tens or 
hundr eds of their member clusters (tFellhauer & Kroupal2002i 
120051) . Furthermore, in more quiescent environments, a sig- 
nificant fraction of clusters may form as bound binaries or 
low-order multiples. Observationally, this is indicated by the 
large number of young binary clust ers observed in the M ag- 
ellanic Clouds (MCs). For instance. iDieball et all fe002l) es- 
timate that about one cluster in eight in the Large Magellanic 
Cloud is a member of a bound group. Most clusters in (can- 
didate) binaries are coeval and younger than 300 Myr, sug- 
gesting that binary clusters generally merge early. Another 
indication that mergers may be common, at least in environ- 
ments such as the MCs, is th at the clusters in these gala xies 
are significantly flattened (e.g.-lKontizasetalJr^S^yr^O^van 
den Bergh 1991J), possibly as a result of rotation that could 
stem from a merger. The typical la ck of significant rotat ion in 
Galactic or M3 1 globular clusters dHan & RvdeiJ!994h does 
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not invalidate the idea that a large fraction of them may also 
be merger products; the more intense tidal field of the Galaxy 
may lead to faster decrease of the rotation. Theoretically, bi- 
nary clusters are predicted to be a common outcome of the 
off-center collision of two molecular clouds (Fujimoto & Ku- 
mai Tl997inB"ekkT'et'lil^ of a spherical 

shell of stars whose f ormation is t riggered by a SN explosion 
in a molecular cloud (lTheisl2002l) . 

In this paper, we consider the possibility that the merg- 
ing clusters each contain an intermediate-mass black hole 
(IMBH; mass M ~ 1O 2_4 M ) and investigate the conse- 
quences for gravitational-wave (GW) astronomy. The obser- 
vational evidence for IMBHs in clusters is only suggestive 
Jvan der Marelll2004l) . However several authors have shown 
how mass segregation leads in sufficiently compact young 
clusters to a phase of runaway collisions among the heaviest 
stars, possibly leading to the formation of an IMBH (Porte - 
gies Zwart & McMi llan l2000HGurkan et alJl2004t Portegies 
Zwart et al. l2004HFreitag et alJ2006blK 

2. FORMATION OF A MERGED CLUSTER 

To form a realistic merged stellar cluster we make two clus- 
ters collide on a parabolic orbit so that the minimum distance 
at which they pass by is d^m if they are considered to be point 
particles at their respective centers of mass. We chose for our 
fiducial simulation d m \ n = 2 pc, corresponding to a relative ve- 
locity at pericenter of 23. 3km s -1 . Given that the clusters are 
set initially on a parabolic orbit, the inevitable loss of energy 
(to, e.g., escaping stars) means that they will always merge to 
form a larger cluster. Each cluster contains AC = 6.3 x 10 4 
particles of 1 M each which are distributed according to 
King models with Wo = 7 (cluster 1) and Wo = 6 (cluster 2). 
Their core radii are R core \ = 0.203 pc and /? C ore2 = 0.293 pc. 
Thus, they are compact enough to have experienced an early 
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FIG. 1. — Evolution of the IMBH binary of the N— body simulation. In 
the upper panel, we show the shrinkage of the semi-major axis as a solid 
line. The initial part of the curve in dashes represents the separation between 
the IMBHs before they have formed a bound binary. In the middle panel 
we display the evolution of the eccentricity, which increases to ~ 0.8. The 
bottom panel gives a an indication of the relevance of GW emission in the 
evolution of the binary in comparison with stellar dynami cal processes : the 
solid line corresponds to the timescale for GWs emission i Peters 1 964 ) and 
the dashed curve indicates the dynamical shrinkage timescale, given by a/d, 
where d is the time derivative of a. One sees that at the end of the simulation 
when«~5x 10~ 5 pc, GW emission is still less efficient than interaction with 
stars by several orders of magnitude. 



core collapse that m ay le ad to IMBH formatio n according to 
iGurkan et all ( 12004 and lFreitag et all J2006bl) who find that 
for clusters with a moderate initial concentration (Wo ~ 6— 7), 
a half-mass radius smaller than 1 - 2 pc is required. We note 
that our results, to be presented ahead, can be rescaled to 
more extended clusters. The central velocity dispersions are 
of CT C o re i = 8.41kms~ and cr core 2 = 8.29kms _1 . Both clusters 
harbor an intermediate-mass black hole (IMBH) with a mass 
of 300M Q , so that the total particle number is AC = 126002. 
The calculations were performed on special-purpose GRAPE- 
6A PCI cards with the direct-summation NBODYA code of 
lAarsethl d2003l) . The particular advantage of this code for 
our study is Aarseth's close-encounter regularization scheme, 
which is free of any softening. The peak performance of these 
cards is of 130 Gflop, roughly equivalent to 100 single PCs, 
which makes possible simulations of stellar clusters with a 
realistic particle number. 

3. EVOLUTION OF THE IMBH BINARY 

In order to have a global view of the evolution of the system, 
in Fig.Q we present the dynamical evolution of the IMBH bi- 
nary (hereafter "binary"). The length and time units used are 
/?nb = 1 pc and ?nb = 4. 16 x 10 4 yr. In Fig. 10 we show both the 
semi-major axis of the binary and its orbital period versus its 
eccentricity. Whilst the direct summation simulations provide 
us with a highly accurate description of the orbital parame- 
ters of the binary as it approaches the LISA band, the required 
computational effort to follow the evolution up to the very fi- 
nal stage of coalescence is unjustified. Therefore in Fig.© we 
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FIG. 2. — Inspiral of an IMBH binary followed in the eccentricity-semi- 
major axis plane. The jagged solid line (in red in the color online version) is 
the result of the iV— body simulation. The dash -dotted lines (in magenta) are 
results of the analytical model (equations [JU for two sets of parameters. In 
one case, we use Kq = 0.5 and start at an early stage, a = 10 _3 pc, e = 0.42. 
The other track is for Kq = 0.1 and starts with the final values of the W-body 
run, a = 5 X 10~ 5 pc, e = 0.8. In both cases the parameters for the stellar 
background -estimated from the /V-body run- are set to a = 23 km s -1 and 
p = 1.6 X 10 5 Mqpc~ 3 . T he smooth so lid (black) lines are trajectories for 
evolution by GW emission i Peters 1964) and the dashed (blue) lines indicate 
the corresponding timescale, ?gw (labels on the left). The dark dashed area, 
at the bottom right, indicates approximately the region of unstable orbits. 
Finally the light (orange) shaded area shows when the n = 2 harmonic of the 
GW signal is in the LISA bandwidth (P olb < 2 X 10 4 s). 



stop the direct-summation calculation after the initial strongly 
fluctuating phase; when the eccentricity achieves a large and 
seemingly steady value of ~ 0.8. We note that an eccentricity 
larger than 0.65 is reached in additional simulations we have 
performed with slightly different initial conditions and that a 
similar result (e ~ 0.8) was found in a simulation of a merger 
of two clusters of 10,000 particles each using a special regu- 
larization method fo r the massive binary (S. Aa rseth, personal 
communication; see Mik kola & Aarsethl2002l for the numeri- 
cal method). In order to determine the binary properties in the 
LISA frequency band, we employ the results of the Af-body 
simulation at t ~ 550 ?nb — 23 Myr and extend it with a ana- 
lytical method down to the moment in which the binary en- 
ters LISA band. The basic idea is to split the evolution of the 
semi-major axis and the eccentricity into two contributions, 
one driven by the dynamical interactions with stars (subscript 
dyn) and the other due to emission of GWs (subscript GW), 
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The "GW" terms ar e as giv e n in i PetersI ( 119641) . Using 
the relationships of lOuinlanl ( 119961) . we have that since 
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FIG. 3. — Characteristic amplitude of the first five harmonics of 
the quadrupolar gravitational radiation emitted during the inspiral of an 
300 Mq + 300 Mq binary. The orbit is evolved according to Eq. \f\ starting 
with the conditions similar to those reached at the end of the fiducial simula- 
tion, i.e., e = 0.8, a = 5 X 10~ 5 pc with an ambient stellar density p c \ and (1-D) 
velocity dispersion cr c i . We assume the sou rce is at a distance D = 1 Gpc. W e 
indicate the noise curve VfSi.(f) for LISA ILarson et at. 2000 ; Lars oii2003> . 
with the Galactic bina ry white dwarf confusion background in dashed line 
IBender & Hils 1 992^and^in^ravVfor the BBO sensitivity jCurve^Cutler 
& Harms 120061) . We label the position of the source at various times be- 
fore plunge. Note that the height of the point for the amplitude above the 
LISAIBBO curve does not represent the SNR (see text). 



In this last expression, G is the gravitational constant, p is the 
density, a is the velocit y dispersion and H is the hardening 
constant jQuinlanll9?6l). We measure H ~ 16, in ag reement 
with previous works (IOuinlanll996tlSesana et alJ2004l) . For 
a hard enough binary, de/d\n{\ /a)|dyn = K(e) and so 



de 
dt 



dyn 



Gp 

■ H—aK(e), 
a 
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with K(e) ~ K e(l - e 1 ) dMerritt & Milosavlievicll2005l and 
references therein). We choose the value of Kq = 0.1 for 
matching purposes with the direct-summation /V-body sim- 
ulation. During the run, the ratio p/a decreases by less than 
30 % at the radius of influence of the binary so we assume it is 
constant. The evolution is dominated by GW emission during 
the last 10 8 yrs (the evolution from the moment the binary is 
bound takes ~ 160 Myr) and enters the LISA bandwidth with 
e ~ 0.07. If we emolov the fiducial value of Kn = 0.5 (Merritt 
& Milosavlievic 120051) but start earlier on the same curve (at 
a lower eccentricity), we find a evolution which fits the initial 
part of the /V-body run satisfactorily but produces too large 
an eccentricity at later times. 

4. DETECTION BY LISA AND THE BBO 

In the late phase of its inspiral, a binary may become a de- 
tectable source of GWs. The characteristic amplitude of the 
gravitational radiation from a source emitting at frequency 
/ is h c = (TrD)~ 1 (2E/f) 1 ' 2 , where D is the distance to the 
source, E is the power emitted and f the time derivative 
of the frequency dFinn & ThorndEOOOl) . With this defini- 
tion, the signal-to-noise ratio (SNR) of an event is obtained, 



assuming ideal signal processing, by the integral (SNR) 2 = 
If'[h 2 c (D/fS h (f)]d(lnf), where /, and f 2 are the initial and 
final frequencies of the source during the observation and 
Sh (f) is the instr umental noise of the detector at frequency 
f llPhinnevl200llBarack & Cutlerl2004l) . 

In Fig. |5] we follow the signal emitted by the IMBHs of the 
fiducial simulation during their GW-driven inspiral. We plot 
the five lowest harmonics of the uuadrunolar emission (Peters 
& Mathews 119631) . In this figure, we assume a distance of 
1 Gpc. Both for LISA and the BBO, only the n = 2 harmonic is 
detectable, during the last few years of inspiral. However, the 
small residual excentricity induces a difference in the phase 
evolution of the « = 2 signal compared to a circular inspiral. 
If the source is followed from a time f mlg before merger until 
merger, the accumlated phase shift is 

17/12 / M _ \ 25/36 
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1 yr 
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where e w -4 Hz is the eccentricity when the n = 2 sig- 
nal has reached a frequency of 10~ 4 Hz and M. z = (1 + 
7)(MiMi ) 3 / 5 (M, +Mor'/ 5 is t he redshifted chiro mass (Cutler 
& Harms 120061) " For our fiducial case (with e I0 -4 Hz = 0.07), 
the eccentricity should be detectable (Atp e > 2ir) if observa- 
tions span at least the last 3-4 years before merger. A mea- 
sured non-zero eccentricity is an important constraint on the 
formation process and stellar environment of the binary. The 
BBO could detect a 300M Q + 300M Q inspiral to at least a 
redshift of 17 but the residual eccentricity will not be mea- 
surable (fnug < 1 day) and significantly more massive IMBHs 
cannot be seen because their GWs are redshifted to too low 
frequencies. 

We now proc eed to estimate rough ly the detection rate of 
binary IMBHs. iFregeau et al.l J2006I) have recently consid- 
ered the detection of binary IMBHs by LISA, assuming that, 
generic ally, two IMBHs form in any cluster undergoing a col- 
lisional run away, as suggested b y the stellar dynamical sim- 
ulations of iGurkan et al.l J2006I) . Here we assume that the 
runaway process leads to the formation of only one IMBH; 
indeed it is not clear whether the two collision-grown very 
massive stars (VMSs, with m* 3> 100M Q ) can avoid merg- 
ing with each other before they become IMBHs. Compared 
to the scenario considered by these authors, ours involve one 
more step, namely, the merger of a IMBH-hosting cluster with 
another one. Schematically, the computation of the detection 
rate for both scenarios can be written 
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(5) 



D m ax ls the distance at which a coalescence event can be de- 
tected, n cmst is the rate of formation per unit volume of clus- 
ters massive enough for (potential) IMBH formation. Rough 
estimates for these terms are Z) max sa 1 — 3 Gpc and n c i ust > 
5 x 10" 10 Mpc" 3 yr 1 ( IFregeau et alJl2006l) . P hin is the prob- 
ability that a cluster will host a (coalescing) IMBH binary. 
Only this term differs between the two scenarios. In our case, 
it can be written as 

p — p p^ p^ 

'bin - 'merg^raMMBH 

where P ra is the probability that a cluster evolves to the run- 
away phase, Pimbh is the probability that the runaway leads 
to IMBH formation and P merg is the probability tha t the cl us- 
ter merges with another cluster. For lFregeau et alJ i2006l) . it 
is simply Pbin = ^ra^MBH- Eq. |5]is an oversimplification be- 
cause the various factors are not independent of each other; 
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for instance, Z) max is larger for more massive IMBHs, which 
probably introduces a dependence on M c \. Also the cosmo- 
logical evolu tion of the cluster for mation rate should be taken 
into account. feregeau et al.l (12006ft present a much more rigor- 
ous computation but we can use the above equations to derive 
our rate estimate from their results, T rj P meig P m Tpregeau with 
r Freg eau« 40- 50(P ra /0. l)yr l . 

We take the same optimistic view as Fregeau et al. that 
VMSs always form IMBHs, Pimbh = 1. The extra P a factor 
indicates that both clusters must contain an IMBH. The con- 
ditions for runaway are relatively well understood (Giirkan 
et al. 12004 iPortegies Zwart et al.ll2004 iFreitag et al.ll2006bt 
iGiirkan et al.ll2006l) but the value of P m is highly uncertain 
because these conditions apply to the initial cluster proper- 
ties. Fregeau et al. choose P la = 0.1 as an illustrative value. 
We note that if all clusters are born with a concentration as 
high as a Wo = 8 King model and the mass and half-mass radii 
of observed local globular clu sters are similar to th eir initial 
values, P m may be about 0.5 llFreitag et al.ll2006cl) . On the 
other hand, iBaumgardt et al.l (120041) have argued that Galac- 
tic globular clusters containing an IMBH may be amongst 
the least dense clusters (for their mass) as a result of the 
strong IMBH-oowered gravothermal expansion (see also Fre- 
itag et al. l2006al) . The relatively small number of observed 
low-density clusters would therefore suggest a smaller P rd but 
the galactic tidal field may prevent the expansion of the clus- 
ter. 

The factor Pmerg is also uncertain. As mentioned in the in- 
troduction, observational and theoretical points support the 
possibility of a large fraction of clusters being born in bi- 
naries, suggesting Pmerg = 0.1 — 1. Furthermore, if the sce- 
nario of Fellhauer & Kroupa holds, each UCDG would re- 
quire ~ 10-100 cluster-cluster mergers. UCDGs may be as 
numerous as normal galaxies, at least in galaxy c lusters (e.g., 
iHilker & Mieske1l200l iDrinkwater et al.ll2005l) . Hence the 
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number of clusters that have been incorporated into UCDGs 
may be of the same order of magnitude as the number of those 
that have survived as isolated clusters (possibly after a binary 
merger). This indicates that the rate estimate for the "UCDG 
channel" should be similar (given the considerable uncertain- 
ties) to that of the "binary-merger channel". Finally, we find 
that the detection rate estimate for LISA is of 

r = 4-5(P ra /0.1) 2 P meig yr- 1 . 

Because of its superior sensibility at frequencies higher 
than ~ 10" 3 Hz, the BBO should be sensitive to IMBH binary 
coalescences from redshifts to at least ~ 15. A high merger 
rate could be envisaged as a potential confusion foreground 
for this miss ion. Fortunately, followi ng similar arguments as 
discussed in ICutler & Harms! d2006h regarding the contribu- 
tion of neutron star mergers, IMBH binaries should be easy 
to subtract out (C. Cutler personal communication). Thus the 
BBO will open a new window of possibilities for our under- 
standing of astrophysical scenarios leading to coalescence of 
IMBH binaries. 

ACKNOWLEDGMENTS 

We thank Leor Barack, Curt Cutler and Jonathan Gair for 
enlightening discussions and Cole Miller for helping us im- 
proving the manuscript. We acknowledge the Astronomisches 
Rechen-Institut for the computing resources on the GRACE 
cluster of the Volkswagen Foundation, SFB439. The work of 
PAS has been supported in the framework of the Third Level 
Agreement between the DFG and the IAC (Instituto de As- 
troffsica de Canarias). The work of MF is funded through the 
PPARC rolling grant at the Institute of Astronomy (IoA) in 
Cambridge. PAS is indebted with the IoA for inviting him the 
month of February 2006, where MF and he started working 
on this project. 



Kontizas E., Kontizas M., Sedmak G., Smareglia R., 1989, AJ, 98, 590 
Kontizas E., Kontizas M., Sedmak G., Smareglia R., Dapergolas A., 1990, 
AJ, 100, 425 

Larson S. L., 2003, Online sensitivity curve generator. 

http://www.srl.caltech.edu/ shane/sensitivity/ 
Larson S. L., Hiscock W. A., Herrings R. W., 2000, Physical Review D, 62, 

062001 

Merritt D., Milosavljevic M., 2005, Living Reviews in Relativity, 8, 8 
Mikkola, S. & Aarseth, S. 2002, Celestial Mechanics and Dynamical 

Astronomy, 84, 343 
Peters P. C, 1964, Physical Review, 136, 1224 
Peters P. C, Mathews J., 1963, Physical Review, 131, 435 
Phinney E. S., 2002, LISA Science Requirements, available on-line at 

www.its.caltech.edu/~esp/lisa/ LISTwgl.req-pr.pdf 
Portegies Zwart S. F, Baumgardt H., Hut P., Makino J., McMillan S. L. W., 

2004, Nat, 428, 724 
Portegies Zwart S. F, McMillan S. L. W., 2000, ApJ Lett., 528, L17 
Quinlan G. D., 1996, New Astronomy, 1, 255 
Sesana A., Haardt F„ Madau P., Volonteri M., 2004, ApJ, 61 1, 623 
Theis C, 2002, Ap&SS, 281, 97 
van den Bergh S., 1991, ApJ, 369, 1 

van der Marel, R. P. 2004, in Coevolution of Black Holes and Galaxies, from 
the Carnegie Observatories Centennial Symposia., ed. L. Ho (Cambridge 
University Press), 37 

Whitmore B. C, Zhang Q., Leitherer C, Fall S. M., Schweizer F„ Miller 
B. W., 1999, AJ, 118, 1551 

Zhang Q., Fall S. M., 1999, ApJ Lett., 527, L81 



